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We explore how a recently developed analytical black-hole binary spacetime can be extended 
using numerical simulations to go beyond the slow-inspiral phase. The analytic spacetime solves 
the Einstein field equations approximately, with the approximation error becoming progressively 
smaller the more separated the binary. To continue the spacetime beyond the slow-inspiral phase, 
we need to transition. Such a transition was previously explored at smaller separations. Here, we 
perform this transition at a separation of D = 20M (large enough that the analytical metric is 
expected to be accurate), and evolve for six orbits. We find that small constraint violations can 
have large dynamical effects, but these can be removed by using a constraint-damping system like 
the conformal covariant formulation of the Z4 system. We Hnd agreement between the subsequent 
numerical spacetime and the predictions of post-Newtonian theory for the waveform and inspiral 
rate that is within the post-Newtonian truncation error. 
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I. INTRODUCTION 

The field of numerical relativity (NR) has progressed 
at a remarkable rate since the breakthroughs of 2005 dl¬ 
ls], when it hrst became possible to simulate the late in¬ 
spiral, plunge, merger, and ringdown of black-hole bi¬ 
naries (BHBs). Recently, Lousto and Healy |3| com¬ 
pleted a long-term 50-orbit processing BHB simulation 
using the moving punctures approach, and Szilagyi et 
al. [5] completed the longest BHB simulation to date: 
the last 176 orbits for a nonspinning, intermediate-mass- 
ratio {mijm 2 = 1/7) BHB using the generalized har¬ 
monic approach. This is a remarkable achievement, but 
the scaling of the inspiral time with the initial separa¬ 
tion T ^ means that evolving a binary through the 
long inspiral is prohibitively expensive, even for highly 
efficient codes. Such a simulation becomes even more ex¬ 
pensive when one is interested in performing long-term 
dynamical evolutions of the relativistic magnetohydrody¬ 
namics (MHD) of circumbinary disks around inspiraling 
supermassive BHBs (SMBHBs). This is because the cir¬ 
cumbinary gas can exhibit significant secular variations 
on the time scale of hundreds to thousands of binary or¬ 
bits. 

In order to make these long-term simulations possible, 
our group developed a complementary approach to treat 
dynamical BHB spacetimes. In a series of papers EHini, 
we used an analytic spacetime that is an approximate 
solution to the Einstein field equations in the inspiral 
regime to describe the evolution of the accretion disks 
surrounding the binary and each of the individual BHs. 

Our initial approach |^, was one in which relativistic 


effects were present but relatively small. In the situation 
when gravity is weak [r^/r = GM/{rc^) I] and mo¬ 
tions are slow [(u/c)^ ^ I], the post-Newtonian (PN) ap¬ 
proximation gives a very good description of spacetime. 
One can then simply construct a PN metric which takes 
energy loss from the binary into account, accurately mod¬ 
eling both the mass loss and inspiral of the binary [Ill- 
Using a spacetime accurate to 2.5PN order [i.e., includ¬ 
ing terms up to ~ but describing the binary 

orbital evolution to 3.5PN, we demonstrated that cir¬ 
cumbinary disks can track the inspiral of a SMBHB un¬ 
til the binary practically reaches the relativistic merger 
regime [5|. The shortcoming of this approach was that 
the PN metric was not valid very close to the BHs, and 
consequently, we excised any material that fell within 1.5 
binary separations. This prevented us from studying the 
dynamics of the gas all the way down through the hori¬ 
zons of each BH. 

In a more recent paper [8|, we extended the metric 
to cover the full BHB spacetime up to the rapid plunge 
state. We did this by extending the framework estab¬ 
lished in Refs. [nfflu for constructing a spacetime metric 
valid for initial data, i.e., a metric accurate for all spa¬ 
tial points but in a very small time interval, to develop 
a metric valid for arbitrary times. In this approach, the 
near zone (NZ), i.e., a zone well outside the two BHs, 
but less than a gravitational wavelength from the binary, 
is still described using a PN expansion. In the far zone 
(FZ), i.e., farther than one wavelength from the binary, 
the metric is described by a post-Minkowskian (PM) ex¬ 
pansion. Finally, near each BH, i.e., in the inner zone 
(IZ), the metric is described using a perturbed Kerr (here 
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Schwarzschild) BH. The metrics covering the different 
zones are smoothly stitched together using asymptotic 
expansions and transition functions in their overlapping 
regions of validity. 

This approach allows us to follow inspiraling SMB- 
HBs over hundreds to thousands of binary orbits, the 
timescale on which gas accumulates, without having to 
solve the Einstein equations numerically. The numerical 
advantage here is that the numerical time step is lim¬ 
ited by fluid characteristic speeds, rather than the much 
faster speed of light (this advantage is diminished if we 
want to evolve gas right near the horizons). 

Here, we explore the possibility of using a hybrid ap¬ 
proach; i.e., use the analytic metric for the long inspiral 
down to separations where our global spacetime metric is 
still valid, and then transition to a full numerical simula¬ 
tion using the analytic spacetime as initial data (the use 
of PN techniques to generate initial data for BHBs was 
first developed in Refs. [HHiH]). The way to do this is 
to convert our approximate spacetime prescription into 
suitable initial data for 3 -I- 1 NR evolutions, evolve the 
data forward in time, and compare the orbital evolution, 
test particle trajectories, and gravitational radiation out¬ 
put with our approximate solution. 

Perhaps more well known is the complimentary ap¬ 
proach of combining PN and other analytical techniques 
with numerical waveforms to generate highly accurate 
hybrid waveforms. Many authors have explored this and, 
we refer the reader to Refs. 0 [IMl and references 
therein. 

The main motivation of this paper is to develop 
techniques to smoothly transition from an analytically 
evolved spacetime to a numerically evolved one. By 
smooth, here we mean that all families of geodesics pass¬ 
ing through the transition region will have continuous 
second derivatives. Here we are considering test particle 
trajectories as stand-ins for fluid trajectories. In partic¬ 
ular, if there is a jump in the second derivative of the 
fluid, we can expect a quasiequilibrium fluid configura¬ 
tion to shock and therefore require a reequilibration that 
may take longer than the inspiral time. Of course, bulk 
binary dynamics are important too. Therefore, we want 
the physics of the inspiral (rate, orbital frequency) to be 
as unaffected by the transition as possible. 

We note that the use of PN techniques to generate 
consistent initial data (i.e., data with the correct radi¬ 
ation content) provides the final ingredient proposed in 
the “Lazarus approach” to generating waveforms 
The proposal there was to transition from PN to NR tech¬ 
niques and then from NR to perturbative techniques. 

When using the global analytic metric as initial data, 
the resulting initial data are essentially equivalent (there 
is only a small difference in the NZ/FZ transition func¬ 
tion and the two metric prescriptions coincide at t = 0) to 
the initial data proposed in Johnson-McDaniel et al. [14] 
and first evolved in Reifenberger and Tichy |33|. Reifen- 
berger and Tichy compared evolutions of Bowen-York 
data [33] to several different analytic initial data con¬ 


structions, including Johnson-McDaniel et al.. Our work 
here extends upon the work of Reifenberger and Tichy 
in several ways, (i) We use the full 4-dimensional metric 
of [S] to compare the dynamics of the numerically evolved 
metric with the analytic one, (ii) we evolve binaries with 
separations large enough that the PN metric and binary 
dynamics are expected to be accurate, and (iii) we find 
techniques to ameliorate the inaccuracies associated with 
evolving these data that were discovered by Reifenberger 
and Tichy. These inaccuracies arise both from constraint 
violations, due to the fact that our global metric solves 
the Einstein field equations only approximately, and from 
inaccuracy in the PN orbital angular momentum and in¬ 
spiral rate, which can lead to eccentricity in the numerical 
binary evolution. 

For the current work, we start the full numerical sim¬ 
ulations when the binary is separated hy D = 20M and 
evolve for six orbits. As shown in Ref. [35], where the 
authors there explored numerical simulations of Bowen- 
York data at separations ranging from D = lOOM to 
D = 20M, there is good agreement between the predic¬ 
tions of PN theory and numerical simulations at D = 
2QM. Additionally, simulations starting at D = 20M 
down to merger are possible with our current codes, as 
demonstrated in Ref. [4] (such a simulation would require 
approximately 1 10® CPU hours on an AMD Opteron 
machine). 

This paper is organized as follows: In Sec.]^ we review 
how the analytic BHB inspiral metric is constructed, as 
well as how it is used to generate 3-1-1 initial data. In 
Sec. jlll] we describe the techniques we used to numeri¬ 
cally evolve the spacetime metric. In Sec. IV we provide 
details on how the simulations were performed and the 
key outcomes of these simulations. In Sec.[Vj we compare 
the results from the numerical simulation at separations 
~ 20M with the predictions of PN theory. Finally, in 
Sec. IVII we discuss our results both in terms of the accu¬ 
racy of the binary dynamics (e.g., inspiral rate and or¬ 
bital frequency) and in terms of gravitational waveform 
generation. 

Throughout this paper, we use the geometric unit sys¬ 
tem, where G = c = 1, with the useful conversion factor 
IMq = 1.477 km = 4.926 x 10"® s. 


II. ANALYTIC BHB INSPIRAL METRIC 

In this paper, we restrict our analysis to non-spinning 
BHs in quasicircular orbits. In this context, it is useful 
to provide a brief review here of our approximate solu¬ 
tion to the Einstein field equations of a BHB spacetime 
in the inspiral regime [8]. The inclusion of spins, both 
aligned [7] and unaligned in this spacetime framework 
will be the subject of future studies. 

This framework was first introduced in Refs. [I2l - H4] as 
initial data for BHB evolutions, and was generalized in 
Ref. [5] to be a full BHB spacetime. In this framework, 
the spacetime is constructed by asymptotically match- 
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TABLE II. The orders of the various approximations used in 
the global metric. The virial theorem implies that mjr is 
taken to be 0{v^ jc?), and the O symbol denotes the highest- 
order term included in the expansion. Here gu and gtt refer 
to the diagonal components of the metric; the rest are the off- 
diagonal components. Note that the version of the first-order 
metric used in Refs. differed from this table in that gu 
was only 0{vfc)'^ there; here, we use resummation techniques 
on both the first- and second-order metrics (in Refs. EHH 
resummation was only used for the second-order metric). 



First order 

Second order 

IZ multipole 

£ = 2 static 

1 = 2, 1 = 2> static 

NZ gu 

0{vyc‘^) 

0{v‘ jc') 

NZ gu 

0{v^/c^) 

Otv^jc^) 

NZ gu 

0(g?lc^) 


NZ g^j 

0 

Otv^jc^) 

FZ gu 

0{v‘^lc^) 

OfgFW) 

FZ gu 

Otv^jc^) 

jc^) 

FZ gu 

Oiv^jc^) 

Oiv^jc^) 

FZ gij 

Otv^jc^) 

Otv^jc^) 


FIG. 1. Schematic diagram of the zones. BHl and BH2 are 
denoted by solid black dots, where the orbital separation is 
ri 2 . The BZs are denoted with gray shells, the outer one rep¬ 
resenting the FZ/NZ BZ and the two inner ones representing 
the NZ/IZ BZs (see also Table|^. The IZ, NZ and FZ are also 
shown in the figure. 


TABLE 1. Location of the inner, near, and far zones, as well 
as the buffer zones joining them. Here rin and rout are the 
approximate inner and outer boundaries of a given zone, mt is 
the mass of BH i, n is the coordinate distance to BH i, ri 2 is 
the binary separation, and A is the wavelength of gravitational 
radiation emitted by the binary. 


Zone 

rin 

?^out 

Region 

IZ BHl (n) 

0 

< ri2 


IZ BH2 (r2) 

0 

< ri2 


NZ {ca) 

» rriA 

< A 


FZ (r) 

IZ-NZ BZ 

> ri2 

oo 

mA<.rA<. ri2 

NZ-FZ BZ 



ri2 <C r <C A 


ing metrics in three different zones characterizing three 
different spacetime regions of validity for different ana¬ 
lytic metrics: (i) a far zone (FZ) where the spacetime 
can be described by a two-body perturbed flat spacetime 
with outgoing gravitational radiation and where retar¬ 
dation effects are fully accounted for; (ii) a near zone 
(NZ) which is less than one GW length from the center 
of mass of the binary (but not too close to each BH) that 
is described by a PN metric (this includes retardation 
effects at a perturbative level and binding interactions 
between the two BHs); and (iii) inner zones (IZs) that 
are described by perturbed Schwarzschild (or Kerr) BHs. 
The full spacetime is then constructed by smoothly tran¬ 


sitioning from zone to zone in the so-called buffer zones 
(BZs). A schematic diagram of these zones and a ta¬ 
ble describing where the zone boundaries are located are 
provided in Fig. and Table |T] (these were also presented 

in Refs. mm- 

In the sections below, we will refer to these initial data 
as the second-order analytical metric. It is constructed 
by asymptotically matching a 2.5PN metric in the NZ 
[the matching is only done for terms up through 0(v/c)'^] 
to a Schwarzschild metric with quadrupole (and its time 
derivatives) and octupole tidal deformations in the IZ. 
As explained in detail in Ref. [5], this matching is ap¬ 
proximate in the sense that it does not lead to a formal 
second-order asymptotic matching in all metric compo¬ 
nents for all times. However, as demonstrated there, it 
does lead to a significant improvement against a lower- 
order analytic metric, the first-order analytical metric, 
which is constructed by asymptotically matching a IPN 
NZ metric [only terms of order 0[v/c)^ are matched] into 
a Schwarzschild metric with quadrupole tidal deforma¬ 
tions. The matching for this first-order metric is exact. 
The metric in the FZ is constructed from the PM ex¬ 
pansion over a flat spacetime with source multipolar de¬ 
composition, where the source multipoles are expanded 
in the PN approximation up to 2.5PN. Note that in the 
PM formalism, the PN metric in the NZ and the multipo¬ 
lar metric in the FZ are formally asymptotically matched 
up to 2.5PN in the NZ-FZ BZ. The precise orders used 
for the calculation of the metric pieces composing these 
analytical metrics are given in Table |TT[ 


III. TECHNIQUES 

We evolved the BHB initial data using the LazEv |3S] 
implementation of the moving puncture approach mm 
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with the conformal function W = ^Jx = exp(—2(?!)) sug¬ 
gested by Ref. [37] and the Z4 [38114(1] and BSSN [411143] 
evolution systems. Here we use the conformal covariant 
Z4 (CCZ4) implementation of Ref. [ 10 ] ■ Note that the 
same technique has been recently applied to the evolution 
of binary neutron stars [S] |45]. For the CCZ4 system, 
we again used the conformal factor W. We used centered 
eighth-order finite differencing for all spatial derivatives, 
a fourth-order Runge-Kutta time integrator, and both 
fifth- and seventh-order Kreiss-Oliger dissipation [46]. 

Our code uses the EinsteinToolkit [TTIHO] / Cac¬ 
tus m / Carpet [OUIOO] infrastructure. The Carpet 
mesh refinement driver provides a “moving boxes” style 
of mesh refinement. In this approach, refined grids of 
fixed size are arranged about the coordinate centers of 
both holes. The Carpet code then moves these fine 
grids about the computational domain by following the 
trajectories of the two BHs. 

We use AHFinderDirect [S3] to locate apparent 
horizons. We also use the Antenna code [33] to calcu¬ 
late the Weyl scalar 

We measure the distance between the two BHs using 
the simple proper distance or SPD. The SPD is the proper 
distance, on a given spatial slice, between the two BH 
apparent horizons as measured along the coordinate line 
joining the two centers. As such, it is gauge dependent, 
but still gives reasonable results (see Ref. [3S] for more 
details). 

To obtain initial data, we use eighth-order finite differ¬ 
encing of the analytic global metric to obtain the 4-metric 
and all its first derivatives at every point on our simu¬ 
lation grid. The finite differencing of the global metric 
is constructed so that the truncation error is negligible 
compared to the subsequent truncation errors in the full 
numerical simulation (here we used finite difference step 
size of 10“"^, which is 90 times smaller than our smallest 
grid size in any of the numerical simulations discussed 
below). We then reconstruct the spatial 3-metric ■jij and 
extrinsic curvature Kij from the global metric data. Note 
that with the exception of the calculation of the extrin¬ 
sic curvature, we do not use the global metric’s lapse 
and shift. In order to evolve these data, we need to re¬ 
move the singularity at the two BH centers. Unlike in 
the puncture formalism [54], the singularities here are 
true curvature singularities. We stuff [55H57] the BH in¬ 
teriors in order to remove the singularity. Our procedure 
is to replace the singular metric well inside the horizons 
with nonsingular (but constraint violating) data through 
the transformations 


f{r) Itj i¥=j, 
-t /(r) 7« + (1 - /(?’))S, 
A„ ^ /(r) A„ , 

where 


f{r) 


0 r < Tmin 

1 r > ; 

P(p) ^min 'C rmax 


( 1 ) 

( 2 ) 

(3) 


(4) 



FIG. 2. The metric function IF at t = 0 plotted versus x 
showing the effects of stuffing the BH interior and the same 
function at t = 165M (but plotted versus y). The more regu¬ 
lar shape of W near the center of the BH at t = 165M is typ¬ 
ical of moving puncture simulations (note that the puncture 
is offset from the y axis by O.OIM). For reference, the func¬ 
tion W for a Bowen-York puncture simulation (solid curve) 
when the puncture crosses the x axis for the second time is 
included (the plot of the Bowen-York data has been shifted). 
Note how the stuffed W appears to evolve to be very similar 
to the standard trumpet W typical of puncture initial data. 


r is the distance to a BH center, and P{r) is a fifth-order 
polynomial that obeys P(rmin) = R'(rmin) = P''{rmin) = 
0, R(rmax) = 1, P'irmB.x) = P”{rnu^x) = 0, and S is a 
large number. The resulting data are therefore glob¬ 
ally. The parameters rmin, ^’max, and S are chosen such 
that both transitions occur inside the BHs and so that 
W varies smoothly with negligible shoulders in the tran¬ 
sition region and is small at the centers. In Fig. [^ we 
show the profile of the conformal factor W at t = OM 
and t = 165M. The former clearly shows the effects of 
stuffing while the latter shows that the system appears 
to evolve to the standard moving puncture gauge (i.e., 
the conformal function W takes on the usual profile for 
a trumpet slicing just like it does when using puncture 
initial data). 


IV. SIMULATIONS 

The initial data parameters for our BHB simulations 
are given in Table |TTT| To evolve the second-order ana¬ 
lytical data, we used the following grid structure: The 
coarsest grid spanned 0 < x < 3200M, —3200M < y < 
3200M, and 0 < z < 3200M (we used 7r-rotation sym¬ 
metry and z-reflection symmetry). The refinement levels 
were centered on the two BHs with half-widths of 1600, 
800, 440, 220, 110, 55, 25, 10, 5, 2, and 0.75. In the figures 
below we denote the resolution of the coarsest grid by h. 
Our lowest-resolution runs had a coarsest resolution of 
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TABLE III. Initial data parameters, mi and m 2 are the 
masses of the two BHs, D is the orbital separation, 5fIorb 
and 5r are the modifications to the 3.5PN orbital frequency 
and inspiral required to reduce the eccentricity, H is a scale 
factor (see text), and th is the measured horizon radius. 


m\/M 

0.5 

m2/M 

0.5 

D/M 

20.00 

M 5Horb 

7.88515 X 10“® 

Sr 

-1.54103 X lO"'^ 

Train/M 

0.05 

^’max/A/ 

0.25 


800 

th/M 

0.484 




h = ho = 32M. We increased the resolution by suc¬ 
cessive factors of 1.2 for the higher-resolution runs. Our 
standard choice, which we used for all long-term runs 
shown below used a medium resolution of h = ho/\.2. 
The highest-resolution run had h = ho/1.2^. 

To demonstrate the smoothness of the transition from 
analytical to numerical evolution we, evolve a set of test 
particles using both the second-order analytic metric and 
the numerical evolution of the second-order analytic data. 
Note that at t = 0, the 4-metrics associated with the two 
evolution schemes are geometrically identical (i.e., they 
only differ by a coordinate transformation implicit in us¬ 
ing different choices for the lapse and shift at t = 0). 
However, because the evolution schemes are different, the 
two spacetimes will have different effective stress ener¬ 
gies (i.e.. Gab will differ for the two spacetimes) even at 
t = 0. Thus, even if the two spacetimes were initially 
in the same coordinate system, higher-order time deriva¬ 
tives (third and higher) of the geodesics will not agree. 
Thus, we only expect continuity of the force acting on 
the geodesics as we transition from analytical to numer¬ 
ical evolutions of the metric. As shown in Fig. we do 
see a relatively smooth transition at early times with the 
two sets of geodesics initially agreeing quite well and then 
deviating more significantly at later times. Importantly, 
this latter deviation is due to two effects, differences in 
the later time coordinate systems and differences in the 
curvature. Plots of the same geodesics from the CCZ4 
and BSSN evolutions are nearly identical. 

Figure provides evidence that the transition from 
analytical to numerical evolutions is sufficiently smooth 
that no sudden impulses are imparted to timelike 
geodesics. However, we still need to demonstrate that 
the subsequent dynamics are accurate. This will require 
that the constraint violations do not significantly affect 
the dynamics of the binary and that the binary remains 
quasicircular. 


A. Initial Hamiltonian and momentum constraint 
violations 

Evolutions of second-order analytical data with BSSN 
were first performed in Tichy |33) . where, like we do 
here, they looked at mass conservation, inspiral time, 



FIG. 3. A set of timelike geodesics initially equally spaced in 
X and normal to the t = 0 hypersurface (note that one of the 
BHs is centered at a; = lOM at t = 0). Here the proper time of 
each geodesic is plotted as a function of the geodesic’s spatial 
position. The solid (black) curves correspond to geodesics 
evolved on the numerical spacetime using BSSN, the dot- 
dashed (cyan) curve correspond to the geodesics evolved on 
the numerical spacetime using CCZ4, finally the dotted (red) 
curves are for geodesics evolved with the second-order an¬ 
alytical metric. The plot to the right zooms in on a typi¬ 
cal geodesic near the start of the simulation. Note how the 
numerical spacetime geodesics smoothly deviate from their 
analytical spacetime counterparts and there is no noticeable 
difference between the test particle trajectories for the BSSN 
and CCZ4 spacetimes at these early times. The rapid change 
in gauge near the start of the simulation is apparent in the 
smooth change in the geodesic seen at a proper time of about 
T = 3M in the inset. 


anomalous eccentricity, and constraint violations (see 
also Ref. [H])- The conclusion there, as well as here, is 
that residual constraint violations lead to relatively large 
errors in the subsequent dynamics. 

Our expectation is that inaccuracies in the second- 
order analytic metric will decrease as the binary sepa¬ 
ration is increased. To test this assumption, we need a 
measure of the constraint violation that is related to the 
dynamics of the binary. Since we can interpret violations 
of the Hamiltonian constraint as an unphysical matter 
field on our spacetime, a natural measure of the degree 
of violation is the total amount of unphysical matter com¬ 
pared to the total amount of physical mass (in this case, 
the total amount of physical mass is « IM). One sub¬ 
tlety we have to contend with is that both positive and 
negative mass densities are dynamically important, and 
there is no reason to expect their respective contribu¬ 
tions to the total error will cancel. Thus, we consider 
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mass 



FIG. 4. The unphysical mass of the binary versus binary 
separation D for the second-order analytical data (i.e., initial 
data) as measured using the integrals munphys. and rriabs. • The 
horizontal axis is in units of D/M, where D is the separation 
of the binary. The larger mass is mabs.. Note that while a 
significant amount of unphysical matter is present, it is spread 
out such that only a fraction of it is absorbed by the BHs (see 
Fig.[5f. Note also that there is a near equal amount of positive 
and negative mass (which is why munphys. is about a factor 
of 20 smaller than mabs.). Rather than plotting munphys., 
we plot —munphys., slncc the net unphysical mass is actually 
negative. 


two measures of the unphysical mass given by 


munphys. - j 

[ IdyD/dv, 

(5) 

""-bs. = ^ j 

[ \‘H\y/jdv, 

(6) 


where H is the Hamiltonian constraint violation and the 
integral is over the cube centered on the origin with side 
length 400M excluding the interiors of the horizons. The 
former gives the net unphysical mass in the spacetime, 
while the absolute value in the latter ensures that positive 
and negative mass densities do not cancel. 

We show both of these measures of unphysical mass 
versus separation in Fig. Note that there is a near 
equal amount of positive and negative mass (which is why 
Wunphys. IS about a factor of 20 smaller than mabs.)- The 
magnitude of munphys. decreases with binary separation 
as roughly {D/M)~^-^, while mabs. decreases at the rate 
of {D. The masses were calculated at t = 0 for 
three resolutions. In all cases, the truncation errors for 
the highest resolution corresponded to an uncertainty in 
the second or higher significant digit in the mass. 

At a separation oi D = 20M, we find that munphys. = 
O.OOIM, while mabs. = 0.033M (Bowen-York data for a 
D = 20M binary solved using the TwoPunctures [58] 
code with 40^ collocation points gives | mabs. I < 7 x 
Note that mabs. increases much more rapidly 
than the power-law prediction with decreasing separation 
for D < 25M, but munphys. is only 30% larger than the 
power-law prediction at D = 20M. It is important to 
note that the BHs will not absorb this much mass, as 


their cross sections are quite small. 

An important result from Fig. is that while the un¬ 
physical mass tends to zero at infinite binary separation, 
for practical purposes, it is never small in the regime 
where we would use NR evolutions. Thus, we need a way 
of removing the unphysical mass from the system. 

We also examine how the quantity of unphysical mat¬ 
ter (mabs.) depends on the locations of the inner/near 
buffer zones. For our runs, we used the transition pa¬ 
rameters of Ref. [14] , which were optimized for a separa¬ 
tion ot D = IQM. Using these parameters, we find that 
mabs. is 0.033M. By optimizing the parameters to reduce 
mabs., we reduce this by only 4% to 0.032M. Interest¬ 
ingly, while mabs. is reduced, the constraint violations 
are more concentrated near the two BHs, thus allowing 
for more absorption of constraint-violating matter by the 
BHs. Importantly, the constraint violations cannot be 
significantly reduced by moving the locations of the zone 
boundaries, since they are fixed analytic functions of the 
masses and separation delimiting overlapping regions of 
validity for different metric approximations. 

It is important to determine not only how much un¬ 
physical matter is present, but also where it is located. 
To this end, we plot the Hamiltonian constraint viola¬ 
tions on the equatorial plane for BHBs at separations 
oi D = 20M and D — lOOM using both the standard 
second-order analytical metric described in Sec. |TT| and 
the first-order version. The main difference between the 
two is described in Table im The Hamiltonian constraint 
violations show a clear improvement as we switch from 
the first-order to the second-order analytical metric, in¬ 
dicating that it is the low PN order which dominates 
the error. Perhaps unexpectedly, even at a separation of 
D = lOOM, the first-order metric has mabs. = 0.02M, 
while the second-order metric has mabs. = O.OOIM. Ex¬ 
amining Fig. we see that the constraint violations are 
concentrated in extended clouds well outside the horizons 
in the buffer zone between the inner and near zones. Most 
importantly, the first-order metric shows a strong shell 
of high constraint violation surrounding the two BHs. 
The second-order metric, on the other hand, has a lower- 
amplitude, more diffuse cloud of constraint violation that 
is less likely to be absorbed. 


B. Effects of the Hamiltonian and momentum 
constraint violations 

In this section, we examine how the numerical evolu¬ 
tion scheme can compound or mitigate issues associated 
with initial constraint violations. To this end, we evolved 
the second-order analytical data using both the BSSN 
formulation pTI - HH] and the constraint-damping CCZ4 
approach. 

Our initial explorations of the dynamics of the second- 
order analytical data were based on the BSSN and CCZ4 
systems. As shown in this section, we find that the CCZ4 
is uniformly better than BSSN in evolving data with non- 
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X-Axis (M) X-Axis (M) 



FIG. 5. The Hamiltonian constraint violation in the vicinity of the two BHs in the binary for binary separations oi D = 20M 
(top) and D = lOOM (bottom). First-order metrics are shown to the left, and second-order metrics to the right. The very 
small white circles at the centers of the BHs are the horizons. 


trivial constraint violations. 

One of the most important differences between the 
BSSN and CCZ4 evolutions is in the horizon mass con¬ 
servation. As shown in Fig. the mass conservation of 
the BHs was relatively poor for BSSN and substantially 
better for CCZ4. In the figure, we see that the BSSN run 
showed an initial increase in mass of 10“^M, followed by 
a mass loss of similar magnitude. While a change of 2 
parts in 1000 may seem small, the effect of this mass 
change on the orbital trajectory is quite large. 

To determine the cause of the lack of conservation of 
the (apparent) horizon mass, we compare the time deriva¬ 
tive of the horizon mass (iMn/dt = dMi^ 2 /dt {Mi = M 2 
by symmetry) with the average value of the constraint on 
the horizons Hh and the flux of constraint violation into 
the horizon Ch (since the spacetime around the two hori¬ 
zons is identical by symmetry, we only plot the constraint 
violation averaged over one of the BHs). We define Hh 
and Ch as 

_ §Hy^dAdB 
§i/ddAdB ’ 

Ch = ~ f C^iiiy/adAdB , ( 8 ) 



FIG. 6. The individual apparent horizon masses for a BHB 
with initial separation of 20M for both GGZ4 (black) and 
BSSN (red) evolutions of the second-order analytical data. 
The BSSN curves show very large oscillations, while the GGZ4 
curves show a much smaller linear growth. The large mass os¬ 
cillations in the BSSN run are due to absorption of constraint 
violations. Note how the effects changing the evolution sys¬ 
tem from BSSN to CGZ4 are much larger than truncation 
error effects. 
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FIG. 7. The Hamiltonian constraint averaged over the hori¬ 
zon as a function of time. This linear log plot shows the 
absolute value of the constraint violation for the BSSN and 
CCZ4 evolutions of the second-order analytic metric, as well 
as the evolution of Bowen-York data (BY) using CCZ4. The 
BSSN, CCZ4, and BY data points are colored according to 
the sign of the constraint violation. The lower plot shows the 
early-time behavior. 


where H is the Hamiltonian constraint violation, C® is 
the momentum constraint violation, rii is the unit (out¬ 
ward) normal to the horizon, y/adAdB is the proper area 
element on the horizon, and the integrals are performed 
over the surface of the horizon. 

In Fig. we show the constraint violation averaged 
over the individual horizons for both BSSN and CCZ4. 
A large positive violation is observed at early times for 
BSSN, which is followed by a negative constraint viola¬ 
tion. This links well with the initial increase in hori¬ 
zon mass for BSSN, which is followed by a later-time de¬ 
crease. The CCZ4 constraints are a factor of 100 smaller 
and do not appear to be correlated with the CCZ4 hori¬ 
zon mass. Because of these results, all of our long-term 
simulations used CCZ4. 

Aside from an overall positive (constant) numerical 
factor, plots of Hh, Ch, and dAdn/dt are nearly iden¬ 
tical for BSSN (see Fig.[^. This means that all three are 
strongly correlating (i.e., dMn/dt (x Ch oc 'Hh)- This 
provides a compelling argument that it is the constraint 
violations that cause the horizon masses to fluctuate. On 
the other hand, for CCZ4, there is no compelling correla¬ 
tion between dMn /dt and the (much smaller) constraint 
violations (see Fig. [^. 

Finally, we examine how the constraint violations in 
the bulk of the simulation domain behave with time. As 



FIG. 8. The Hamiltonian constraint averaged over the hori¬ 
zon, the flux of the momentum constraint violations through 
the horizon, and the time derivative of the horizon mass for 
the BSSN simulation. The Hamiltonian and change of rate of 
the horizon mass have been multiplied by constant positive 
factors. Note the near-perfect correlation of horizon mass 
change and constraint violation on the horizon. In the hgure, 
solid curves correspond to a resolution of /io/1.2, while dotted 
curves correspond to a resolution of ho/1.2^. 



FIG. 9. The Hamiltonian constraint averaged over the hori¬ 
zon and the time derivative of the horizon mass for the GGZ4 
simulation. Note how unlike in Fig.[^ these two are not corre¬ 
lated. Here both the constraint violation and m appear to be 
converging to a very small value but from opposite directions. 


shown in Fig. 10 the norms of the constraint viola¬ 
tions for CCZ4 and BSSN evolutions behave quite differ¬ 
ently (here we restrict the norm to the volume inside 
a ball of radius SOM and outside the two horizons so 
that the norm is dominated by constraint violations rel¬ 
atively close to the binary). The CCZ4 constraints fall 
to a much lower level (about a factor of 1000 smaller for 
the Hamiltonian and a factor of 100 for the momentum) 
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FIG. 10. Time evolution of the norms of the constraint 
violations for BSSN and CCZ4 evolutions. Note how the 
BSSN Hamiltonian constraint remains relatively high while 
the CCZ4 constraints quickly fall to about 10“^. Also shown 
are the norms of the constraint violations for an evolution 
of Bowen-York data with CCZ4. The Bowen-York data leads 
to constraint violations that are on average a factor of 2 below 
constraint violations for the new data. 


than BSSN. For comparison purposes, we also performed 
equivalent evolutions with standard Bowen-York initial 
data [34] . 

As shown in Fig. El we see clear convergence of the 
norm of the momentum constraint violation to zero 
for CCZ4, while the Hamiltonian constraint violation, 
though small, seems to bottom out at about 10“®. The 
Hamiltonian constraint for an equivalent Bowen-York run 
bottoms out at roughly half this value. For this conver¬ 
gence check, we only ran the highest-resolution run for 
350M due to computational costs. 

The amount of unphysical mass that the BHs can ab¬ 
sorb depends not only on the amount of unphysical mat¬ 
ter, but also on the dynamics of the unphysical matter. 
For BSSN, the constraint violations largely stay in place 
(and can therefore be accreted) due to the presence of a 
zero-speed constraint mode in BSSN [39] , while for CCZ4 
the constraints quickly leave the vicinity of the BHs. The 
very different behavior of BSSN and CCZ4 is shown in 
Fig.jg 

The overall efficacy of using CCZ4 to drive the con¬ 
straint violations to zero can be measured by examining 
in detail how well the horizon masses are conserved. As 
shown in Figs. |13| and |14[ there is a relatively strong linear 
trend in the mass that, while converging to a small value, 
is substantially larger than the Bowen-York result. Here 
we also see a significant advantage to using higher-order 
dissipation. Note that even with the highest-resolution 
runs, the horizon mass increase is an order of magnitude 
larger than for Bowen-York data evolved with CCZ4. 
Since the Bowen-York data were evolved with the same 
evolution system and grid structure, it appears that there 


FIG. 11. The norm of Hamiltonian constraint violations 
and the norm of the Euclidean magnitude of the momen¬ 
tum constraint violations. The momentum constraints have 
been rescaled by a factor of (ho/h)'^, where h is the base res¬ 
olution of a particular run and ho is the base resolution of 
the coarsest simulation. Note that the Hamiltonian viola¬ 
tions have not been rescaled. The spikes occurring roughly 
every 350M are due to the high-frequency initial gauge wave 
reflections off the mesh refinement boundaries. Here too, the 
norms have been restricted to the interior of a sphere of 
radius SOM about the origin and outside the two horizons 


are peculiarities associated with the analytic initial data 
driving the mass increase (note, as seen in Fig. [9 absorp¬ 
tion of constraint violation seems not to be the cause 
of this mass increase). One possibility which we have 
not explored in detail here is that the methods used to 
stujf the BHs may be affecting the mass conservation. 
The other candidate would be residual constraint viola¬ 
tions. Even for our highest-resolution run, the average 
constraint violation on the horizon surface at late times 
was 50% larger than for Bowen-York (the global con¬ 
straint violations were, on average, an order of magnitude 
larger, as shown in Fig. 10). As observed in Ref. |9j, dif¬ 
ferences in accuracy of the spacetime at this level will 
likely not be important for MHD simulations. 


C. Eccentricity 

As shown in Fig. the transition from analytical to 
numerical evolutions is relatively smooth, and, as shown 
above, evolutions with CCZ4 drive the constraint viola¬ 
tions down to acceptable levels (i.e., within a factor of 
10 of the levels obtained by evolving the constraint sat¬ 
isfying Bowen-York data). The last step required for a 
successful continuation of the evolution is to ensure that 
the binary remains quasicircular (the PN inspiral used 
to generate the data is quasicircular). To accomplish 
this, we applied the eccentricity reduction procedure of 
Ref. |S1] to our data (we found that we needed to set 
MSQorh = 7.88515 x 10"® and 6r = -1.54103 x lO""*). 





























































10 



X-Axis (M) X-Axis (M) 

FIG. 12. The Hamiltonian constraint on the xy plane near the two BHs at t = 200M for the BSSN (left) and CCZ4 (right) 
evolutions of the second-order analytical data starting at a separation of = 20M. The scale is logarithmic and goes from 
10~^ to 10“®. The absolute value of the constraint violations for the BSSN simulation match closely the violations on the 
initial slice, while the CCZ4 violations are three orders of magnitude smaller. In each plot the interiors of the BHs have been 
masked out. The constraint violations at t = 0 are given in the top right plot in Fig. 



FIG. 13. The mass of the individual horizons versus time for 
the CGZ4 simulations using the standard fifth-order dissipa¬ 
tion, seventh-order dissipation, and fifth-order dissipation of 
Bowen-York data. The linear trend in the mass, while con¬ 
verging to a small value, is substantially larger than that for 
Bowen-York. 



FIG. 14. The rate of horizon mass increase versus time. Here 
D5 indicates that fifth-order dissipation was used, while D7 
indicates that seventh-order was used. 


After three iterations, we were left with a residual eccen¬ 
tricity of e = 0.002, which was small enough for this test 
(see. Fig. 15). In Fig. 15 we show the SPD versus time 
for both CCZ4 and BSSN evolutions of the eccentricity- 
reduced data, we also show a CCZ4 evolution of the orig¬ 
inal data. From the figure, the nonphysical dynamics 
(overall increase in radius) of the BSSN evolution is ap¬ 
parent. Note that we implement the eccentricity reduc¬ 
tion by changing the initial orbital inspiral rate f and 
orbital frequency florb used in the PN equations of mo¬ 
tion. Changes to the inner-zone and far-zone metrics are 
automatically handled by the matching procedure. 


The eccentricity reduction here is complicated by the 


fact that the amount of constraint-violating fields ab¬ 
sorbed by the BHs changes as the trajectories are mod¬ 
ified. This in turn, leads to a more complicated depen¬ 
dence of the eccentricity on the orbital parameters than 
is seen for constraint-satisfying data. 


D. Waveform 

The waveforms presented below are relatively short 
(due to the expense of running the simulation to merger, 
which would take about six months on 80 Opteron cores). 
We will be comparing the numerical waveforms to PN 
waveforms. For these “short” runs, the dominant error 






























11 



FIG. 15. The SPD versus time as calculated using a CCZ4 
evolution of the new data both before and after the eccentric¬ 
ity reduction procedure. For reference, a BSSN evolution of 
the low-eccentricity data is also given. Constraint violation 
drives the eccentricity for this latter case. 



FIG. 16. The {£ = 2, m = 2) mode of the waveform (ri/) 4 ) 
extracted at r = 800M and r — 1600M, a linear extrapolation 
of these to r = oo, and an extrapolation to r = oo using 
perturbative techniques. Both the standard extrapolation in 
r and the perturbative approach give very similar waveforms. 


is due to finite extraction radius. In Fig. |16[ we show a 
“late” segment of the waveform extracted at r = 800M, 
r = 1600M, and a linear extrapolation (in I = 1/r) to 
r = oo from these two waveforms, as well as the ex¬ 
trapolation to r = oo using the perturbative approach of 
Refs. [nnilSI]- As expected, the dominant errors due to 
finite radius are phase errors. Finally, in Fig. |17[ we show 
the waveform (post-initial burst) extracted at r = 800M 
for the resolutions h = /io/1.2^ and h = ho/1.2^. We 
also show the extrapolation of these waveforms using the 
techniques of Refs. [iniiii]- As can be seen, the domi¬ 
nant error in the waveform is the phase error due to finite 
extraction radius. 


V. COMPARISON TO PN 



To gauge the accuracy of our transition from analyti¬ 
cal to numerical evolutions, we compare the subsequent 
dynamics of the binary with the predictions of PN. 

In Fig. [T^ we show the SPD versus time and PN sep¬ 
aration versus time. Since the SPD at t = 0 is larger 
than 20M we translate the SPD vertically. Note that 
the SPD is not expected to be equal to the PN separa¬ 
tion. The SPD includes effects due to the nonflatness 
of the spatial metric and measures how distant the two 
horizons are from each other, while the PN separation 
extends from the center of one BH to the other and 
the proper separation corresponding to this would not 
be finite. While it is interesting that the numerical SPD 
matched the PN separation reasonably well, these are not 
gauge-invariant quantities. We also show the SPD for an 
equivalent Bowen-York simulation (evolved with BSSN) 
first reported in Ref. [3S]. The Bowen-York run also had 
the eccentricity reduction procedure applied. 


FIG. 17. The (£ = 2, m = 2) mode of the waveform {r%jj 4 ,) 
extracted at r = 800M at two resolutions. Both the raw 
waveform and the extrapolation to infinity are shown. The 
dominant error is the extrapolation error which manifests it¬ 
self predominately as a phase error. 


To have a more gauge-invariant measure of the accu¬ 
racy of the evolution, we compare the waveform (as ex¬ 
tracted at 1600M) with the 3.5PN prediction for quasi¬ 
circular orbits [52] (similar to what was done in Refs. [55] 
and [24]). All waveforms are shown in Fig. 

When extracting at r = 1600M, we get very good 
agreement between the raw {£ = 2, m = 2) mode of tpi 
and the extrapolation to infinity using the techniques of 
Ref. [60] . Note that the numerical waveform prior to the 
burst of radiation is purely a function of the initial data. 
The initial data used the 3.5PN equations of motion; thus 
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FIG. 18. The SPD and PN separation versus time as calcu¬ 
lated using a CCZ4 evolution of the new data and an older 
BSSN evolution of Bowen-York data. The SPD is always 
larger than the coordinate (and hence PN) separations. We 
shift the SPDs downward by 3.25M so that they agree with 
initial PN separation at t = 0. 


the agreement in the frequency at early times with 3.5PN 
is expected. On the other hand, the PM metric in the 
FZ contains terms up to 2.5PN order, which naturally 
leads to a lower-order approximation for the wave ampli¬ 
tude, since it depends not only on the orbital parameters, 
but also on the metric perturbation order. After the ini¬ 
tial data burst, the waveform becomes noisier, but the 
agreement with 3.5PN is still quite good. The numerical 
waveform amplitude, however, seems to be closer to the 
average of 1.5PN and 3.5PN. 

One important note is that the PN waveform given in 
the initial data is slightly out of phase with the resulting 
numerical waveform, as shown in Fig. That is to say, 
after translating the waveform in time by r* (the tor¬ 
toise coordinate of the extraction observer), the PN and 
NR waveforms agree quite well for the part of the wave¬ 
form after the initial burst has hit the extraction sphere 
(in the plot, this would be from t = 0 to t = 2400M). 
However, prior to this burst arriving at the observer (in 
the plot, prior to t = 0), the PN and NR waveforms are 
out of phase by 0.255 radians. This initial part of the NR 
waveform is produced by the far-zone metric in the initial 
data, while the latter part of the waveform is produced 
by subsequent fully nonlinear binary dynamics. This will 
have repercussions if one wants to smoothly attach a PN 
waveform to the numerical waveform. It is important to 
note that other than a translation by the tortoise coor¬ 
dinate r* corresponding to the extraction radius, the NR 
and PN waveforms have not been translated. 

The phase error itself can be explained by how we con¬ 
struct the metric in the far zone. In the far zone, the 
metric at some point at a (coordinate) distance r from 
the origin depends on the dynamics of the binary at a re¬ 
tarded time given by the light propagation time from the 




FIG. 19. The frequency and magnitude of the (^ = 2, m = 2) 
mode of ril )4 as measured at r = 1600M in the full numerical 
simulation and various PN predictions. Here we use either the 
1.5PN or 3.5PN expressions for r and either the full 3.5PN 
expression for /122 (as a function of r and lo) of Faye et al. |62| . 
or we truncate to IPN order. We use the 3PN expression for 
u in all cases. 


binary to that point. We use the expression = t — r, 
which is the flat-space retarded time. A more accurate 
expression would include the mass of the spacetime. For a 
Schwarzschild BH, in harmonic coordinates, the retarded 
time would be 


j.Sch. 

^ret 


= t — 


(r + M) + 2M log 


r + M 
2M 



(9) 


where M is the mass of the spacetime. We thus find that 
for a given waveform frequency w, using the flat-space re¬ 
tarded time will introduce a phase error of approximately 


6(j) = oj 


M + 2M log 


f r + M 



( 10 ) 


Since the binary’s orbital period here is MHorb ~ 0.01, 
and the (£ = 2, m = 2) mode of the waveform has twice 
this frequency, we expect a phase error introduced by 
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FIG. 20. The numerical waveform from the new data, a nu¬ 
merical waveform from equivalent Bowen-York data, and the 
PN waveform. Here the numerical waveforms are shifted by 
r* (the tortoise coordinate at the extraction radius), and the 
PN waveform is unshifted. The phase agreement is good af¬ 
ter t = 0 but breaks down prior to the initial data pulse (for 
the new data) despite the fact that the NR frequency is in 
close agreement with PN for the whole waveform. The jump 
in phase between the early and late part of the waveform is 
likely due to the use of the flat-space retarded time in con¬ 
structing the early waveform. 

the flat space retarded time of « 0.287 radians, which 
is reasonably close to our measured phase error of 0.255 
radians. 

One final note concerns the amplitude of the initial 
data pulse in the waveform. As seen in Fig. |21[ the ini¬ 
tial pulse of radiation is suppressed relative to equivalent 
Bowen-York data. At r = 1600M, the suppression is 
roughly a factor of 2, while at the r = 400M extraction 
radius, the suppression is closer to a factor of 3. This 
is mostly due to the fact that we have initial data which 
model the astrophysical BHB system better and therefore 
possess less spurious radiation content when compared to 
the conformally flat BY initial data. Here, because of the 
resolution of the grid at r = 1600M, the high-frequency 
gauge pulse (near t = —500M) is completely dissipated 
away. The high-frequency components of the initial data 
pulse are similarly suppressed. 


VI. DISCUSSION 

In order to perform accurate GRMHD simulations of 
gas accreting onto a BHB, including the minidisks around 
each BH, we need a spacetime that is accurate for the en¬ 
tire lifetime of the binary, i.e., from the slow inspiral at 
extremely large separation all the way to merger and re¬ 
laxation of the remnant BH. Our goal, therefore, is to 
produce a four-dimensional metric that is accurate out¬ 
side the horizons at all times, of sufficient smoothness 
that timelike geodesics vary smoothly (i.e., are (7°°) ex- 



t/M 

FIG. 21. A plot of the real part of the {£ = 2,m = 2) mode 
of rxlj 4 (shifted in time by —r) for Bowen-York and the hy¬ 
brid data here (denoted by PN). The waveforms extracted at 
r = 400M show a high-frequency pnlse near t — — lOOM dne 
to an nnresolved gauge wave. This high-frequency pulse is 
dissipated away and not visible in the waveform extracted at 
r = 1600M. 


cept at a single transition time where they are C^, and 
the subsequent binary dynamics should match PN pre¬ 
dictions to a high degree in the vicinity of this transition 
time. To do this, we extended an analytic metric that 
is accurate when the binary’s separation is U ^ M by 
continuing the evolution using fully nonlinear numerical 
techniques for closer separations. 

The main questions we addressed here concerned how 
we can accurately transition from using an analytically 
evolved spacetime metric to a fully nonlinear numerically 
evolved metric that describes the binary during the rapid 
plunge and merger. At the transition, we used the an¬ 
alytical spacetime metric to construct initial data for a 
subsequent numerical evolution (as was previously done 
in Reifenberger and Tichy |33j ). Our work builds upon 
Reifenberger and Tichy in two main ways. We start from 
an analytic spacetime that can be extended arbitrarily 
far into the past, and we can thus compare dynamics of 
particles pre- and post-transition. We also perform the 
transition at a binary separation of D ~ 20M, where 
the binary’s dynamics are still well described by PN the¬ 
ory and errors introduced in the gas dynamics by the 
approximate metric are washed out by MHD turbulence 
(see Ref. [5] for an analysis of MHD evolutions on this 
analytical background for various separations). 

In order for the transition from an analytical evolution 
to a numerical one to be smooth enough, the binary’s 
orbital dynamics could not change significantly as a re¬ 
sult of the transition. The binary’s dynamics in the fully 
nonlinear numerical simulation had two main sources of 
error. First, constraint violations led to rapid unphysical 
oscillations in the orbital separation. Second, small errors 
in the PN expressions for the orbital angular momentum 
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and inspiral rate led to eccentricity in the binary. We 
were able to ameliorate the first source of error by evolv¬ 
ing with the constraint-damping CCZ4 [30] formulation 
of the Einstein equations, which causes constraint vio¬ 
lations to rapidly propagate away from the BHs, signifi¬ 
cantly reducing unphysical binary dynamics. In addition, 
by adding small changes to the initial inspiral rate and or¬ 
bital frequency, we significantly reduced the eccentricity 
of the numerical binary using the techniques of Ref. |59j . 


We subsequently found that the NR evolution leads to 
the expected gravitational waveform, orbital frequency, 
and binary inspiral rate (to within the truncation error 
of the simulation). The remaining error we found is a 
phase error in the early part of the waveform. This phase 
error is about 0.255 radians. We ascribe this error to our 
use of the flat-space retarded time in the far-zone. By 
not including effects due to the mass of the spacetime we 
generate phase errors of the order of 0.287 radians in the 
waveform. This error can itself be ameliorated by using 
the Schwarzschild retarded time when constructing the 
far-zone metric, which is something we will explore in an 
upcoming paper. 
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